Reconciling Semiclassical and Bohmian Mechanics: 
V. Wavepacket Dynamics 
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In previous articles [J. Chem. Phys. 121 4501 (2004), J. Chem. Phys. 124 034115 (2006), 
J. Chem. Phys. 124 034116 (2006), J. Phys. Chem. A 111 10400 (2007)] a bipolar counter- 
propagating wave decomposition, 9 — \&+ + \I r _, was presented for stationary states ^ of the one- 
dimensional Schrodinger equation, such that the components ^± approach their semiclassical WKB 
analogs in the large action limit. The corresponding bipolar quantum trajectories are classical-like 
and well-behaved, even when has many nodes, or is wildly oscillatory. In this paper, the method is 
generalized for time-dependent wavepacket dynamics applications, and applied to several benchmark 
problems, including multisurface systems with nonadiabatic coupling. 



I. INTRODUCTION 

Trajectory interpretations of quantum mechanics have 
existed since the beginning of the quantum theory — 
even predating the Schrodinger equation itself. In- 
deed, one such approach survives today in the form of 
the Jeffrey- Wentzel-Kramcrs-Brillouin (JWKB) approx- 
imation, or more generally, semiclassical mechanics i 2 ^ 3 - 
In this approach, a time-evolving quantum wavepacket 
is treated as a statistical ensemble of classical tra- 
jectories that "carry" approximate quantum informa- 
tion, i.e. complex amplitudes. In the early 1950s, 
hearkening back to the earlier pioneers, D. Bohm and 
coworkers developed a conceptually similar trajectory 
interpretation of the exact quantum theory4>^> 6 ' intro- 
ducing the so-called "quantum potential" to guide the 
resultant quantum trajectory dynamics. In the in- 
tervening years, quantum trajectory methods (QTM)s 
have been used "analytically" — to provide insight into 
solved time-dependent quantum wavepacket propagation 
problemsii&SJii — and more recently, as a "synthetic" 
tool — to actually solve the time-dependent Schrodinger 
equation (TDSE) itself^^^^^^^^ Note that in 
this paper, "QTM" refers to the original quantum trajec- 
tory method based on standard Bohmian mechanics, as 
developed by Wyatt and coworkers jii as well as the vari- 
ous offshoots and approximations that have developed in 
the intervening years. 

The standard Bohmian formulation uses an amplitude- 
phase decomposition of the wavefunction, ip, which in one 
dimension (ID) takes the form 



ip(x,t) = R{x,t)e lS{x ' t)/h . 



(1) 



This representation has been called "unipolar,"— be- 
cause the field functions, R{x,t) and S(x,t), are single- 
valued at all positions, x, and times, t. The field func- 
tions, and the resultant quantum trajectories, are gener- 
ally smooth and classical-like — provided the true poten- 
tial, V(x), is slowly- varying, and^{x,t) exhibits no inter- 
ference. However, interference introduces non-classical- 
like oscillations in R(x,t) and S(x,t), which in turn lead 
to severe numerical difficulties for QTM calculations — 



collectively referred to as "the node problem." 11 ' 16 De- 
spite substantial progress ) 20 ' 21 ' 22 ' 23 the node problem 
continues to be the most formidable roadblock imped- 
ing the progress of QTM's as a general and robust tool 
for exact quantum scattering applications. 

As a promising remedy to the node problem, this paper 
is the fifth in a serie a 24 ' 25 ' 26 ' 27 exploring the use of bipolar 
decompositions of the wavefunction, i.e. 



■0 = ^+ +V>-> 



(2) 



such that the quantum trajectories associated with the 
bipolar wavefunction components, ij)±(x,t), are well- 
behaved and classical-like, even when the unipolar quan- 
tum trajectories associated with ip(x,t) itself are not — 
e.g., when ijj(x,t) exhibits interference. In fact, by gen- 
eralizing Bohmian mechanics for multipolar decomposi- 
tions such as Eq. ([2]) above, it becomes possible to achieve 
classical correspondence — i.e., trajectories and field func- 
tions that approach their classical counterparts in the 
classical limit of large mass, energy, and/or action — 
which is not in general possible for any unipolar Bohmian 
treatment. 

Most QTM papers found in the literature concern 
themselves with time-dependent localized wavepacket dy- 
namics, as indeed, is also true of this fifth paper in 
the bipolar series. However, this paper represents a 
stark departure from the previous four, all of which per- 
tain to stationary solutions of the the time- independent 
Schrodinger equation — albeit obtained in a pseudo-time- 
dependent manner using delocalized counter-propagating 
wave components. For ID stationary states, the two 
ip±(x,t) components correspond to approximate semi- 
classical analogs^* 3 - of which there are always two at 
every classically-allowed point in space and time — thus 
justifying a bipolar QTM treatment in this context. If 
there are classical turning points demarcating classically 
allowed and forbidden regions, these are known a priori, 
and do not change over time. 

In contrast, the semiclassical treatment of localized 
wavepacket dynamics, even in ID, is considerably more 
complicated than for stationary states — rendering it far 
more difficult to arrive at a suitable QTM analog. In par- 
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ticular, the time-evolving semiclassical field functions can 
become multivalued over regions of x and t in a highly 
non-trivial way, due to caustics that form, move, and 
then disappear over timei The appropriate number of 
semiclassical components can become one, two, or even 
three or more, and in any case, varies nontrivially over x 
and t. The formation of caustics, in turn, is due to the 
crossing of neighboring classical trajectories belonging to 
a semiclassical wavepacket ensemble. In stark contrast 
are quantum trajectories, which, for a single component, 
may never cross — implying single- valuedness for all x and 
t, and evidently greatly complicating the task of achiev- 
ing classical correspondence. 

Let us imagine, for instance, that at the initial time, 
to, ip(x,t ) is taken to be of the unipolar Eq. ([T]) form. 
Though initially single- valued, ip(x,t) may, at later t, 
become multivalued over certain regions of x, under a 
semiclassical treatment. However, since an ensemble of 
quantum trajectories does not develop caustics, how then 
is a QTM to become similarly multivalued? This funda- 
mental difficulty persists even under the usual allowance 
that quantum and semiclassical wavefunctions are equiv- 
alent only to 0(K). In fact, it raises the following, 
purely semiclassical conundrum, which to the author's 
best knowledge, has not been previously addressed: if 
i/j(x, t) can become multivalued at later t, why should it 
never be considered so at t = to? One practical reason 
is that the 0(K) uniqueness of the semiclassical repre- 
sentation would be compromised, making it unclear how 
to proceed. This answer is not satisfying in any formal 
sense; yet indeed, semiclassical theory treats the same 
Gaussian wavefunction as single-valued if interpreted in 
a wavepacket context, or double-valued if regarded as a 
harmonic oscillator ground state. 

With regard to QTM wavepacket dynamics, the mul- 
tivalued problem described above can in principle be 
addressed in a variety of ways. One simple strategy 
would be to actually use classical trajectories for the 
dynamics — not approximately, as in semiclassical theo- 
ries, but rather, with exact propagation of the (complex- 
valued) quantum amplitudes, R{x, t). This approach can 
be regarded as an arbitrary Lagrangian-Eulerian (ALE) 
method ) 11 ' 21 ' 22 modified to allow for the formation of 
caustics and multivalued fields. Note that probability 
is no longer conserved along trajectories — and indeed, 
must approach zero at the caustics, in order to avoid in- 
finite probability density. In practice, this condition itself 
causes numerical instabilities, and in any case would be 
unfeasible to implement for large systems^ 

Alternatively — and seemingly contrary to semiclassical 
behavior — one might opt to adhere strictly to a globally 
bipolar decomposition of the Eq. @ form, in order to 
facilitate comparison with the bipolar stationary state 
theories of papers I-IV in the series i 24 ' 25 ' 26 ' 27 Even in 
this relatively restricted context, however, the "correct" 
bipolar wavepacket generalization is not necessarily ob- 
vious. In the previous work, for instance, both of the 
ip+ and field functions are symmetric for station- 



ary bound states^ whereas only S is symmetric [i.e., 
S-(x) — — S+(x)] for the more general case of stationary 
scattering states! 25 ' 26 For bipolar wavepacket dynamics, 
neither field function provides us with simplifying sym- 
metry; moreover, we evidently have no direct recourse 
to semiclassical mechanics, which was previously relied 
upon as a guide. 

One way to state the problem is as follows: a complete 
specification of the bipolar wavepacket dynamics gener- 
ally requires four independent real- valued time-evolution 
equations. Two of these are automatically provided by 
the TDSE, but how are the remaining two to be cho- 
sen? One promising avenue, which we have explored 
considerably* 2 ^ is to adopt the combined flux continu- 
ity condition^ 6 .^^ [Eq. fT7])] as the third equation. In 
fact, if imaginary flux is considered,— the fourth and fi- 
nal equation can also be obtained. Unfortunately, these 
equations do not provide satisfactory results, in that over 
time, the individual ip±(x,t) components themselves de- 
velop interference oscillations and nodes — thus defeat- 
ing the purpose of a bipolar expansion. More flexibility 
can be obtained by dropping the imaginary flux conti- 
nuity condition, but to date, all such efforts have also 
been unsuccessful 2 ^ — either due to ip±i x ,t) interference, 
or other equally unsatisfactory behaviors (Sec. Ill C|) . 

Ultimately, the most successful bipolar wavepacket 
generalization scheme we have considered has also proven 
to be one of the most conceptually straightforward — i.e., 
to expand ip(x,t) as a superposition of stationary state 
solutions, whose + and — components are then used to 
determine the wavepacket '4 ! ±{x,t) via linear superpo- 
sition. The idea is simple, but the theoretical devel- 
opment is somewhat involved. In any event, this ap- 
proach, which will serve as the focus of this paper, leads 
to a node-free, and otherwise remarkably well-behaved 
Eq. (J2) decomposition — and moreover, turns out to sat- 
isfy classical correspondence after all. It also leads to 
time-evolution equations that are practicable for numer- 
ical implementation. Note that bipolar quantum trajec- 
tory formulations for the TDSE have been considered 
previously by other authors^ albeit not in a manner 
designed to solve the interference/node problem. 

The remainder of this paper is organized as follows. 
Sec.|TT]first presents the requisite background of the bipo- 
lar theory for stationary scattering states (Sec. Ill A[) . 
followed by a derivation of the new wavepacket time- 
evolution equations for asymptotically symmetric poten- 
tials (Sec. Ill Bp . Additional properties of the resultant 
ip±(x,t) component wavepacket dynamics are described 
in Sec. Ill CI Generalizations for asymptotically asym- 
metric and multisurface applications are then provided 
in Sees. Ill Dl and III El respectively. Results and discus- 
sion, for benchmark applications of each type described 
above, are then presented in Sec. lIIII Finally, a summary 
and concluding remarks are given in Sec. IIVI 
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II. THEORY 
A. Background 

Consider the two bipolar components, (fi±(x), associ- 
ated with a ID stationary scattering state solution of 
the Schrodinger equation, (f> E (x) — 4> E {x) + 4> E (x), with 
energy, E, and left-incident boundary conditions. The 
solution components, <f> E (x), are exact quantum analogs 
of a type of semiclassical JWKB approximation resulting 
from the "generalized Froman" approachi 2 i 26 i 27 i 30 Note 
that for the remainder of this paper, <j) is used in the 
context of stationary state wavefunctions, whereas "0 is 
reserved for localized wavepacket dynamics. The solu- 
tion components 4>±(x) behave as (right/left) traveling 
plane waves in both asymptotes. As x sweeps through 
the potential interaction region, the solution component 
4> E (x) varies smoothly from incident to transmitted wave, 
whereas 4> E {x) varies smoothly from reflected wave to 
zero. 

In several previous articles , 25 i 26 i 27 i 30 an extremely ac- 
curate, efficient, and robust ID numerical algorithm was 
developed for computing 4>±(x), and thus <f> (x). The 
algorithm is a time-dependent relaxation method, for 
which the initial <\> E = <fr E is a plane wave. Over time, a 
reflected wave, 4>-(x,t), comes into being through inter- 
action region coupling due to the potential energy, and 
eventually, <fi E (x,t) relaxes to the true stationary scat- 
tering solution. 

From Ref. H Eq. (12), the <p±{x,t) time-evolution 
equations are 



dt 



+ : UE-V)<}>l- l ~Vct> E (3) 



where primes denote spatial differentiation, m is the 
mass, p — \ / 2mE is the magnitude of the asymptotic 
momentum, and asymptotically symmetric potentials are 
presumed [V(x) — > as x — > ±oo]. The initial value con- 
ditions are given by 



(t>+°(x) = <f> E {x,t ) 



exp 



ipx 
~~h~ 



iEt 



(4) 



4>-°{x) = (t> E (x,t o )=0, 



where t = to — > — oo is the initial time. At the left 
and right coordinate limits, i.e. x = xl — > — oo and 
x = xr — ► +oo, respectively, the boundary conditions 



4>+(XL,t) 



exp 



ipx L 



iEt 



(5) 



(j) E (x R ,t) = 



In general, Eq. (J3J does not satisfy the TDSE, in that 
(d(j) E /dt) ^ -{i/h)H<t> E for all times t, where H is 
the usual Schrodinger equation Hamiltonian. This is 
consistent with the interpretation of this approach as a 



"revelatory" or "relaxation" method ) 25 i 26 ' 30 but implies 
that Eq. itself cannot be used as a basis for deriv- 
ing wavepacket time-evolution equations, for which the 
TDSE must be satisfied at all t. 



B. Wavepacket time-evolution equations 

We therefore consider the asymptotically large time 
limit, t — > +oo, in which Eq. ([3]) not only satisfies the 
TDSE, but relaxes to the exact stationary solution, so 
that 



dt 



-E(j) E and 



dt 



-E^l 



as t — > +oo. 

(6) 

Substituting Eq. ([6]) into Eq. ([3]) and rearranging yields 
the following time-independent expressions for the spatial 
derivatives of the solution (t — > +oo) <j) E (x): 



lE 1 , * 

6± =± r p< 



>±T T [-)V 



(7) 



Equation ([7]) above is consistent with Ref. HOEq. (12), 
with v = p/m — y / '2E/m. Differentiating Eq. ([7]) with 

respect to x, and then using Eq. |0 to substitute for (fc E 
in the resulting right hand side yields: 



t? v ' ^ • h\p 

Equation [5] can then be used to obtain 



(8) 



(9) 



Substituting Eq. (J6j> for E(f)± above then results in the 
following, new time-evolution equations: 



4 



i - e V 

dt =-n H ^ T 2p- 



(10) 



A subtle shift has occured in the transformation from 
Eq. ([3]) to Eq. (| 10[) . First, the latter manifestly satisfies 
the TDSE at all times t. Second, what constitutes the 
coupling contribution (i.e. the last term) in Eq. Q is 
not equivalent to that of Eq. (TlT)|) . Most strikingly, the 
latter coupling is proportional to V rather than to V 
itself. This implies that the coupling vanishes in both of 
the coordinate asymptotic limits, xl and xr, even when 
V(x) is not taken to be asymptotically symmetric. This 
represents quite an improvement over Eq. ([3]) , which can- 
not be applied to asymmetric potentials because there is 
asymptotic coupling in at least one coordinate limit i 26 ! 30 
Although Eq. (JTUJ) can be used with asymptotically asym- 
metric potentials, there are nontrivial ramifications for 
wavepacket dynamics fSec. HIDj) . The final, and perhaps 
most important, observation that will be made regarding 
the coupling contribution to Eq. (|10p is that it is directly 
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proportional to the JWKB quantity denoting the classi- 
cal limit — i.e., |AV'|, where A = 2nfi/p is the de Broglie 
wavelength.— In particular, the classical limit is obtained 
when |AV'| < p 2 /2m — i.e., (V'/p) — > 0. According to 
Eq. (|10p . the coupling vanishes in the classical limit, in 
which the <fi± themselves approach Schrodinger equation 
solutions. In this manner, classical correspondence is es- 
tablished for the new time-evolution equations. 

All of the above still refers to delocalized left-incident 
stationary scattering states <\> E , at definite energies E. 
For every positive E value, the methodology uniquely 
determines a corresponding <j> E and <f> E . In generalizing 
for (left-incident) localized wavepacket dynamics, an em- 
inently sensible strategy is to decompose the wavepacket 
ip as an orthonormal expansion in the stationary states 
(j) E — an expansion which, in turn, is applied to the <f>± 
components themselves, to uniquely determine ip± and 
the Eq. @ bipolar decomposition. Thus, 

/>oo 

ip(x,t) = / a(E)<f> E (x,t)dE and (11) 



to the following wavepacket time-evolution equations: 



ip±(x,t) = / a(E)4>%(x,t)dE. 
Jo 

In principle, the above equations enable ip±(x,t) to be 
completely determined — provided the initial wavepacket, 
ip°(x) = ip(x,to), is specified, and all of the solution 
4>±(x, i)'s are known a priori. In practice, the lat- 
ter requirement defeats the purpose of doing localized 
wavepacket dynamics — i.e., to avoid explicit calculation 
of the delocalized 4> E states. 

Consequently, we apply the Eq. (jlip expansion to 
both sides of Eq. (p~0|) . in order to directly derive time- 
evolution equations for iJj±(x, t). We require that the ex- 
plicit integrations over dE be tractable, so that <j) E and 
<f> E not appear explicitly in the final results for dip±/dt. 



This in turn requires that the right-hand-side of Eq. (|10D 
exhibit no explicit dependence on E or p, which — due to 
the coupling term — is seen not to be satisfied. To make 
progress, we use the identity, 



ip 

n 



(12) 



obtained from Ref. [2J Eq. (9), or from Ref. |3fJ Eq. (11), 
or by using Eq. ([7]) to add cfif' + <j) E> . Substituting the 
integral of Eq. (|12p with respect to x into Eq. (fTT)|) then 
yields 



dt 

where $±(x) 



H4>i ± y w 
^{x')dx', 



$f) 



,(13) 



apart from a term proportional to 4> E (x — > — oo) that 
vanishes when integrated over E via Eq. (jlip (Dirich- 
let wavepacket boundary conditions, Sec. IIIip . Since 
Eq. (|13[) above has no explicit dependence on E or p, ap- 
plying the Eq. (fTTj) expansion to Eq. (fT5|) leads at once 



dj>± 
dt 

where ^±{x) 



ff^± ± (*+ - *_) ,(14) 
ip±(x')dx' (15) 



Equation (|14| can be directly integrated over time, to 
determine the dynamics of the bipolar wavepacket com- 
ponents, ip±(x, t). In addition to the usual TDSE Hamil- 
tonian contribution, Eq. (|14|) includes a coupling contri- 
bution that is proportional to V' (and independent of 
mass). As in the case of Eq. (TT51) . this implies that the 
asymptotic x — > ±oo coupling vanishes, even for asym- 
metric potentials. Since V(x) also vanishes asymptoti- 
cally, we thus find that ip±(x — > ±oo,f) evolves under 
free-particle propagation. Note that throughout the x 
coordinate range, the ip± coupling terms are equal and 
opposite, so that (dtp/dt) = ~(i/h)Hip. Thus, Eq. 
satisfies the TDSE at all t. 



C. Additional properties 

Having defined a set of wavepacket time-evolution 
equations [Eq. (|14p]. we next consider whether these 
give rise to well-behaved tp±(x,t) components at all x 
and t. In general, a great range of behaviors are pos- 
sible for nonstationary state dynamics in ID, many of 
which are undesirable. We therefore first stipulate that — 
apart from exhibiting interference — ip(x,t) itself be well- 
behaved. By this we mean that ipixit) is normalized 
to unity and well-localized at all times t, consisting 
of a single left-incident wavepacket at the initial time 
t — to, and of well-separated left- and right-moving re- 
flected and transmitted branches, respectively, at the fi- 
nal time t = tt — > +oo. Under these assumptions for 
the ip(x,t) wavepacket dynamics, we define well-behaved 
components, ip±(x,t), as those that satisfy the following 
three conditions: 

• Condition 1: Perfect asymptotic separation at to 



and t 



• Condition 2: Well-localized ip±(x,t) and ~$>±(x,t) 
at all t. 

• Condition 3: Node-free components, ip±{x,t), at 
all t. 

Condition 1. means that the initial left-incident 
wavepacket consists solely of ■0+ — i.e., ^p%(x) — ip°(x) = 
ip(x, to), and (x) = 0. It also means that at the asymp- 
totically large final time tf, ip + (x) — ip+(x,tf) becomes 
the right-moving transmitted branch of tpf (x) — tp(x, tt), 
and (x) = ip- [x,t*) becomes the left- moving reflected 
branch. Condition 2. is straightforward, and absolutely 



essential, e.g., for multidimensional generalizations. Con- 
dition 3. is expected to hold at all t — particularly inter- 
mediate times, where ip(x,t) itself may exhibit substan- 
tial interference and/or nodes. 

For the remainder of this subsection, asymptoti- 
cally symmetric potentials are presumed, as defined in 
Sec. Ill Al At the initial time to, the incident wavepacket 
is localized far to the left of the potential interac- 
tion region, so that V(x) is effectively zero, and the 
Eq. ([5]) plane wave boundary condition accurately de- 
scribes <j>+(x, t) — over the whole asymptotic region where 
|V' (a:)| 2 is significant. Thus, the initial Eq. (jTTJ) station- 
ary state expansion is essentially identical to a Fourier 
expansion — or equivalently, the momentum-space repre- 
sentation, ip° (p) . Through the identification p — \/2mE, 
we find that only the p > states contribute in the 
Eq. (fTTjl expansion for tp+(x); the p < states give rise to 
a left-moving contribution, which is presumed to be zero 
at t = to- Condition 1. thus requires that the negative 
momentum contribution to ip°(x) be vanishingly small, 
i.e. 

f \j°(p)\ 2 dp^0. (16) 

This is a very reasonable requirement to impose on (x) , 
for it implies that the initial wavepacket is completely 
incident upon the scattering potential center. 

At the final time tf, the a(E) expansion coefficients 
from Eq. (jTTJ) are identical to their initial values at 
to. Moreover, the reflected and transmitted wavepacket 
branches are localized far to the left and right, respec- 
tively, of the interaction region, so that once again, the 
4> E {x,t) are effectively plane wavesi 26 ' 30 However, these 
asymptotic plane waves are no longer characterized by 
the standard unit normalization of Eq. (0, but must 
instead be weighted by the independent reflection and 
transmission amplitudes, R(E) and T(E), when deter- 
mining the Fourier components, 4> (p)- Note that since 
4>-(xR,t) — [Eq. ©J, there can be no tp- component 
in the right asymptote, implying that the right-moving 
transmitted branch at tf must consist only of a ip+i x ) 
contribution. Similar arguments can be used to demon- 
strate that the left-moving reflected branch at tf must 
consist only of ipL(x). We thus find that Condition 1. 
above is formally satisfied at both to and tf. 

Regarding Condition 2., here again we make use of 
Eq. (fTo) . If Eq. (fitl)) were not satisfied, then ^(p = 
0,i) would in general be nonzero. Since ^{xR,t) — 
V2nh^jj(p = 0,i) [from Eq. (TT5|) ]. the right-asymptotic 
value of ty(xji,t) would approach a constant (in x), 
nonzero value. Thus, \& (ar, t) would be delocalized — even 
if ip(x,t) itself were localized, as it is initially. Over 
time, moreover, due to the coupling term in Eq. ([15]), 
the tp±(x, t) would themselves eventually become delocal- 
ized, even if ijj{x, t) itself were not — clearly, an untenable 
situation, in violation of Condition 2. Conversely to the 
above scenario, the fact that Eq. (fT^) is true implies that 
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(a) 




position x 

(b) 






momentum p 

FIG. 1: Schematic indicating properties of initial wavepacket, 
ip°: (a) position space; (b) momentum space. The solid line in 

(a) represents p°{x) = \i/j°(x)\ 2 , the initial wavepacket den- 
sity, taken to be a Gaussian centered far to the left of the 
interaction region. The dotted line in (a) represents | (x) | 2 , 
where ^(x) = J* ip(x') dx', as per Eq. [fl5]l. Note that ^(x) 
is not localized, owing to the small constant-valued tail that 
extends towards x = +oo. This is because the Fourier trans- 
form, tp°(p), does not satisfy Eq. (|16[1 perfectly — as seen in 

(b) , a plot of |^°(p)| 2 . 



ty±(x,t), and therefore tp±(x,t), are localized — not only 
at asymptotic times, but at all times. Thus, Condition 
2. is also formally satisfied. 

In practice, Eq. (|16[) is never perfectly satisfied (even in 
the to — > — oo limit), but is only approximately correct, 
to some desired level of numerical accuracy. The true 
^f±(x,t) will have a small delocalized constant- valued 
tail — extending to the right towards x — > oo for the def- 
inite integration convention of Eq. (|15p . as indicated in 
Fig. [TJ but non-vanishing no matter which integration 
limits are adopted. For reasonable ip°(p) distributions 
however, the magnitude of these tails can easily be made 
arbitrarily small (and therefore insignificant) simply by 
shifting to xp°(p — po) for sufficiently large po- 
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Finally, we address the all-important Condition 3. 
Though we can offer no formal proof at present that this 
condition is always satisfied, it must certainly be true 
in the classical limit, in which (V'/p) approaches zero. 
More generally, if Condition 3. is satisfied for the in- 
dividual 0± components — as has been demonstrated for 
a great range and variety of test application a 25 i 26 i 27 i 30 — 
then it is reasonable to expect this property to also be 
preserved under the Eq. (jlip expansion. In any event, 
Condition 3. is verified for each of the test cases consid- 
ered in Sec. 

On the other hand, one nice property of the (f>± under 
Eq. Q that is definitely not extended to the tp± under 
Eq. (I14| is that of combined flux continuity. In other 
words, unlike all previous bipolar formulations for sta- 
tionary state dynamicS ) 26 i 27 ' 30 we find here that 



dp+ dp_ 
dt 



ot 



(17) 



where p± — \ip±\ 2 are component densities, and j± — 
p±(S±' /m) are component fluxes, defined in the standard 
quantum manner via 



ip± (x,t) = R± {x, t) exp [iS± {x, t )/Ti] . 



(18) 



The appearance of the standard flux in the above ex- 
pressions represents a departure from the stationary state 
formalism — for which predetermined classical velocities, 
±w, are used rather than (S'±/m) — and a move back to- 
wards standard Bohmian mechanics On the other 
hand, the fact that the Eq. (|17[) combined flux conti- 
nuity condition is not satisfied might lead one to argue 
that the (S±/m) velocity field is inappropriate here, and 
that some other choice — perhaps some nontrivial gener- 
alization of the stationary state velocities — would lead to 
an Eq. (|17p -tvpe equality. In fact this is incorrect — for it 
can be shown fSec. HII A II) that (p+ + P-) dx is not 
conserved over time, implying the Eq. p7[) inequality re- 
gardless of the particular choice of velocity field. In any 
event, the standard Bohmian velocity field is naturally 
obtained when Eq. (fT4"]) is used to derive time-evolution 
equations for the component densities: 



dp± 
dt 



V 



= -i± ± -r Im [V4 (*+ - *-)] , (is) 



From Eq. (fl9]) above, Eq. (fTT]) is easily ob- 
tained. Note that despite the Eq. (jTTJ) inequality, 

(p° + + p°-) dx = + /) ** = i*zp** - 

1, so that globally over time, the (p + + p_) probability 
is conserved. 



D. Asymptotically asymmetric potentials 

Our next task is to generalize the previous discussion 
for the case of asymptotically asymmetric potentials. To 
be completely general, we allow the left asymptotic value, 



Vl = V(xl), and right asymptotic value, Vr = V(xr), to 
be completely arbitrary — i.e., Vl ^ Vr, and neither Vl 
nor Vr need be zero. In this context, it is straightforward 
to generalize Eq. ([3]) for either of the two asymptotic po- 
tential values, but not both simultaneously. Essentially, 
this is done by adopting V (x) = V° as the effective 
potential used to generate classical trajectories ; 26 ' 2 ^ 30 
where the constant V° is chosen to be either V° = Vl or 
V° = Vr. In either case, for the generalized Eq. ([3]), i.e. 



dt 



= T^<P± +UE-V-V°)cj ) 1 ± 



(V - V )^ 

m 

[with p — y // 2m(E — V )], coupling vanishes in one x 
asymptote, but not the other. 

As explored in previous paperS ) 26 i 30 two natural reme- 
dies for the asymptotic coupling dilemma are considered: 
(1) define a smoothly varying effective potential V eS (x) 
such that V eS (x L /R) — V l /r; (2) define a discontinu- 
ous transition at an intermediate dividing point xr>, so 
that V cS {x) =V L + {Vr - V L )Q(x - x D ), where 6() is 
the (heaviside) step function. With respect to deriving 
wavepacket time-evolution equations as per Sec. Ill El op- 
tion (1) poses severe difficulties, in that it is not clear how 
to recouch the relevant equations^ to avoid explicit de- 
pendence on E and/or p. Option (2) on the other hand, 
is straightforward, as we now demonstrate. 

The key property is that Eq. ([20]) — whether for V° = 
Vl, V° = Vr, or an arbitrary V° value — leads to exactly 
the same wavepacket evolution equations as for V Q = 0, 
i.e., Eq. (fH|) . However, the resultant ip± components are 
^-dependent. Thus, the V° = Vl components ipL±, and 
the V° = Vr components ipR± > constitute distinct bipolar 
decompositions, each satisfying Eqs. ([2]) and (fT4|) . This 
can only be true provided the initial (and final) condi- 
tions are different, i.e. ipL± 7^ ^Pr± an d ^ ipR± : 
which in turn implies that Condition 1. (Sec. Ill Cj) must 
be false. In fact, we still find that ip L+ (x) = ip°(x) and 
tpL-( x ) = 0) t> u t both ip R± (x) 7^ 0. Similarly, at tf, the 
transmitted branch of ipf(x) equals il) R+ {x) [no ip R -{ x ) 
contribution], but when expanded instead in terms of 
i/jl±(x), includes nonzero contributions from both. The 

reflected branch of ipf (x) equals ip^_(x). 

In accord with option (2) above, it is natural to define 
a single bipolar decomposition, tp± , that satisfies all three 
conditions of Sec. Ill CI by "gluing" together the asymp- 
totic solutions at the dividing point, xd, as follows: 

tf)±(x,t) = Q(XD - x)tpL±(x,t) + Q(x - XD)lpR±(x,t) 

Such a procedure is analogous to other dividing sur- 
face methods, commonly used in reactive scattering 
applications i 34 i 35 It is not known at present how to time- 
evolve Eq. (f2"Tj) directly, i.e. without recourse to separate 
calculations for ipL± and ipR± over all x and t. How- 
ever, the latter is straightforward to achieve in practice — 
requiring, in addition to Eq. (|14p. only the specific initial 
conditions ipR±(x), which are derived below. 
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First, in the Eq. (fTTj) expansion, the lower limit of the 
integration must be replaced with E m - m = max (Vl, Vr), 
as reactive scattering does not occur at energies below 
E m - m - For ip L +( x ) ~ r( x )i this expansion is still equiv- 
alent to a Fourier expansion, except that the minimum 
allowed (left) momentum value is 



Pu 







y/2m(V R - V L ) 



if V R < V L ; 
otherwise. 



(22) 



Thus, the upper limit in Eq. (|16|) must be replaced with 
Pmin rather than 0, in order that Condition 2. be satis- 
fied. In any event, both ip°(p), and the a(E) expansion 
coefficients in Eq. fjl 1|) . can be computed explicitly from 
ip (x) via straightforward Fourier transform. 

The next step is to relate the 4 > r± decomposition for 
the stationary state solution, <fi E , to the corresponding 
4>l± decomposition. Applying Eq. ([7]) to (j)^ ± , and rear- 
ranging, we obtain 



=F 



ih 

PR 



(23) 



where pr = y / 2m(E — Vr), For purposes of expanding 
the initial wavepacket tp°(x) = tp L +( x )' we can replace 
4> E in Eq. (|23[) with <fif + - Moreover, the x range of in- 
terest is restricted to the left asymptote, where the <j) E ± 
are plane waves of the Eq. ([5]) form, but with p replaced 
by Pl — \j2m(E — Vl)- Making these substitutions in 
Eq. (23|) leads to 



±[1± P l/pr] </>!+, 



(24) 



from which ipR± ( x ) can be obtained via straightforward 
inverse Fourier transform. In particular, Eqs. (111|) and 
(1241 lead to 



a(E) 



1 



1 ± 



E-V L 
E-Vr 



x exp 



y/2m(E~ V L )x - Et 



(25) 



dE. 



Similar arguments can be used to justify the other initial 
and final conditions for ipL± and ipR±, as discussed earlier 
in this subsection. 

Using the explicit initial value conditions of Eq. (|26|) . 
propagation of the ip R ±(x,t) thus becomes as straight- 
forward as for ipL±(x,t). Once achieved, Eq. (|2ip can 
then be used to construct a bipolar ip±(x,t) decompo- 
sition that satisfies all three conditions of Sec. Ill CI On 
the other hand, from a purely practical standpoint, little 
harm would result if one were to simply use ipL±(x,t) 
throughout x and t. The reason is that, unlike Eq. (f2"TJ)) . 
Eq. (|14p exhibits no coupling in either x asymptote. 
Thus, at both to and tt for instance, ipL±(x,t) [and 
ipR±(x,t)] evolve according to free particle propagation, 
which introduces no nodes or interference. Condition 3. 



is therefore satisfied. In practice, this is far more impor- 
tant than Condition 1., which — as discussed above — is 
not satisfied for the if)^ , decomposition in the transmit- 
ted branch of i/j'(x). 



E. Multisurface generalization 

Like the bipolar stationary state theory, the theory of 
bipolar wavepacket dynamics can also be generalized for 
ID multisurface applications. Let / denote the num- 
ber of electronic states considered. A diabatic-like time- 
independent matrix Schrodinger equation is presumed, of 
the form 



H-<j) h = E4> L 



(26) 



where {4>f, 4> E , ■■■<fif} comprise the vector components 
(associated with each of the / diabatic states) of the nu- 
clear stationary state wavefunction, cj> E , and 



H 



1 2m 



d 2 



(27) 



are the components of the / x / Hamiltonian operator 
matrix, H, with i < f and j < f labeling diabatic states. 

The Vij(x) — Vj y i(x) are the diabatic potential en- 
ergy curves, with the i ^ j case denoting the coupling 
potentials. In order to ensure that intersurface cou- 
pling vanishes in the asymptotic limits (required to ob- 
tain asymptotic scattering waves with correct boundary 
conditions) ! 26 i 27 i 30 W e must have Vi^j(xL) = Vi^j(xR) = 
0. However, the asymptotic values for the diagonal po- 
tentials, Vi y i(x), are allowed to be completely arbitrary, 
and in particular, need not be symmetric. Left and 
right asymptotic values are denoted Vil = Vi^(xL) and 
ViR = Vj,i(xR) , resp ectively. 

As per Sec. IIIDI rather than work with generic effec- 
tive potentials Vf s (x) that smoothly interpolate between 
ViL and ViR values^! we instead choose constant effec- 
tive potentials, V cS (x) — V°, with V° arbitrary for now. 
Note that in general, V° can be chosen to coincide with 
at most one of the ViL and ViR, so that we expect no 
more than one of the 2/ component asymptotes to man- 
ifest perfect asymptotic separation (Condition 1. from 
Sec.ira. 

From Ref. [13, the <j>f — <j>f + + 4>f_ components satisfy 



(28) 



where p = •v/2m (E — V°). By combining Eq. ([26]) with 
Eq. j2Hl), we obtain 



J.E 



T ■ 



E 

(29) 
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the multisurface generalization of Eq. (|7|) . Note that the 
same constant p is used for all components i, and in both 
asymptotes Xl and xr. 

By differentiating Eq. (|29|) with respect to x, and oth- 
erwise applying the procedure described in Sec. Ill Bl we 
obtain time-evolution equations for the stationary state 
components, 



at 



■±[ 



H 



kE -r- 



2p 



E 

3^3 ' 



(30) 



the multisurface generalization of Eq. I|10p . Finally, sub- 
stitution of the integral of Eq. ([28]) (with respect to x) 
into Eq. (fBT)]) . and integration over E 1 via an Eq. (fTTj) -tvpe 
expansion yields the following multisurface wavepacket 
time-evolution equations: 



dip 



'i± 



dt 



i 

where 



£[■ 

4=1 



i? 



,(3i) 



ipj±(x')dx' (32) 



Note that Eq. (|3Tj) satisfies the multisurface TDSE, in 
that 



(33) 



For a given value of V , it is a straightforward mat- 
ter to propagate all of the ipi± wavefunction components 
over time, using Eq. (|31j) in conjunction with initial value 
conditions discussed below. The resultant ipi±(x,t) will 
in general satisfy Conditions 2. and 3. from Sec. Ill CI but 
not Condition 1. Note that the particular ipi = ipi+ + ipi_ 
decompositions obtained depend on the value of V°, even 
though the time evolution equations [Eq. (f3"Tj) ] are V°- 
independent. As described in Sec. Ill Dl the V° depen- 
dence manifests in the initial conditions, ^ ± (x). 

In addition to the conventions and conditions al- 
ready adopted for "well-behaved" wavepacket dynamics 
in Sec. Ill CI let us further presume for the multisurface 
case that the initial wavepacket is left-incident on sur- 
face i = 1. Then, V , ($>i)±( a ') = 0, but ip° ± (x) depends 
on V° via an Eq. (f26]) -tvpe expansion (with ipR± replaced 
with tpi±, Vl replaced with V\l, and Vr replaced with 
V°). One natural choice for V° is V° — V\l itself, lead- 
ing to ipi + (x) — ipi(x) and t/ji_(x) = 0. This choice 
leads to perfect asymptotic separation (i.e. Condition 
1.) for ipi±(x, t) in the left asymptote, but generally not 
for ipi±(x,t) in the right asymptote, nor for any of the 
ip<i>X)±{x, f)'s in either asymptote. 



We again reiterate that from a practical numerical per- 
spective, Condition 1. is not required — i.e. perfectly sen- 
sible results for all ipi±(x,t) may be obtained using the 
V° = Vil choice above, or any other reasonable V° value. 
On the other hand, if one is determined to have perfect 
asymptotic separation for both left and right asymptotes 
for all i, this can also be achieved — via introduction of 
a dividing point xd, and the multisurface generalization 
of Eq. (|2Tj) . In effect, this would require that up to 2/ 
separate calculations be performed, corresponding to all 
of the distinct possibilities for V° — ViL and V° = Vm- 
For each of these calculations, Condition 1. is guaranteed 
for (at least) one tpi component in one asymptote, which 
is then singled out in that particular calculation e.g., 
ip2R+ (x > xd , t) , from the V° = V2R calculation. Finally, 
we note that for the asymptotically symmetric special 
case considered in Ref. [2J, where ViL — ViR = for all i, 
then the single choice V° — leads to perfect separation 
in both asymptotes for all components ipi . 



III. RESULTS 

We have applied the bipolar wavepacket time-evolution 
equations derived in Sec. |TT] to a variety of model ID 
applications. The primary goal is to validate numeri- 
cally that this approach satisfies the three conditions of 
Sec. Ill CI especially Condition 3. Consequently, little at- 
tention is paid here to numerical efficiency, and only the 
simplest algorithms are employed, using Eulerian fixed 
grids with uniform spacing in x and t. No trajectories or 
quantum potentials are computed; these will be consid- 
ered in later papers, that will actually solve the TDSE 
by synthesizing quantum trajectories "on the fly." 

In this paper, Eqs. (fH)) and (f3"Tj) are integrated 
over time using the standard first-order forward Eulcr 
method, with fixed time step size, Eulerian fixed 

grids are used to discretize the spatial coordinate x, with 
uniform spacing Ax, and left and right grid edges, xl 
and xr, respectively. Condition 2. from Sec. Ill CI implies 
that Dirichlet boundary conditions, /(xl) = /(x_r) = 0, 
are employed, where /(x) represents any wavefunction 
component or its spatial integral. The spatial deriva- 
tives implicit in the H contribution to Eqs. ([H]) and 
([3Tj) are evaluted numerically using standard symmet- 
ric (two-sided) second-order finite differenced The spa- 
tial integrations are evaluated using closed Newton-Cotes 
formulas — specifically, the two-point trapezoidal rule for 
the second grid point from the left, and the three-point 
Simpson's rule for all other interior grid points^ 

The initial wavepackets are all taken to be of the stan- 
dard Gaussian form, 



■2l 



1/4 



exp [— 7(x — xq) 2 ] exp 



/ ip x 



(34) 

from which ^°(p) can be determined analytically. In all 
calculations, the parameters 7 and po are chosen such 
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that the Eq. (p~6|) integral is negligibly small, i.e. compa- 
rable to the desired level of numerical accuracy for the 
calculation, which is 10 -6 . Similarly, all of the other 
numerical parameters, A, Ax, xl, xr, and Xo, are con- 
verged to the same level of accuracy. For the model ap- 
plications considered here, typical converged parameter 
values in atomic units are as follows: A sa 0.1; Ax « 0.08; 
X L/R — T35. Unless explicitly stated otherwise, the mass 
is taken to be m = 2000 a.u. Computer animations 
(.wmv file format) for all of the wavepacket dynamics cal- 
culations presented in this paper are available as EPAPS 
supplements^ and by direct request from the author. 



A. Eckart barrier system 

The canonical model scattering system for the asymp- 
totically symmetric special case is the Eckart barrierj 38 i 39 
defined via 

V(x) = Vbsech(ax) 2 , (35) 
and specification of the parameters Vb, a, and m. 



1. proton-like mass 

For the first Eckart application considered here, the 
parameter values are chosen to be Vb = .0024, a = 2.5, 
and m = 2000, respectively, in atomic units. This is sim- 
ilar to what has been called the "Eckart A" system in 
previous papersi 26 i ' 30 In atomic units, the parameters 
describing the initial Gaussian wavepacket of Eq. (13~4")) 
are taken to be 7 = 0.35, xo = —7.0, and p = ~3. 28634 
a.u. [Fig. Ufa)]. The p$ value corresponds to an inci- 
dent kinetic energy of .0027 a.u., which is slightly above 
the barrier energy Vb, so that substantial reflected and 
transmitted branches of ipf(x) are obtained. 

The tp±(x,t) components are propagated using 
Eq. (|14p. and the numerical methods described above, 
to a final time i max = (tt — to) — 11600 a.u., where 
reflected and transmitted branches are found to be well- 
separated in left and right asymptotic regions, respec- 
tively. Fig. [2] shows the resultant wavepacket dynamics 
for tp±(x, t) and ip(x, t) at four representative time slices, 
including t = t^ = and t = tf = t max . Although 
only four time slices are presented, the bipolar decom- 
position has been carefully inspected at all intermediate 
times (Fig. \S§ , to ensure that Fig. [2] captures all of the 
relevant dynamics. 

We find that Condition 1., Condition 2., and above 
all, Condition 3. (from Sec. Ill C[) are indeed well sat- 
isfied at all times. Both bipolar components tp±(x,t) 
are remarkably smooth, localized, non-oscillatory, and 
Gaussian-like at all times — despite the fact that ip(x,t) 
itself displays a substantial amount of interference at in- 
termediate times [dashed curves, Fig. [2Kb) and (c)]. In 



addition to interference, the proton-like mass is suffi- 
ciently small that the wavepacket dynamics also mani- 
fest substantial dispersion, as well. In fact — apart from a 
small "spur" that develops on its left side fSec. IIII A 2[) . 
and the fact that its integrated probability decreases 
over time — the tp+ (x, t) evolution resembles that of a free 
particle Gaussian wavepacket, simultaneously translating 
and dispersing its way through the interaction region, 
and in the process, smoothly transforming from the inci- 
dent wavepacket into the transmitted branch of the final 
wavepacket. 

In contrast, the ip- i x , t) wavepacket dynamics — 
though similarly smooth and well-behaved — exhibit 
somewhat different behavior, due to the different ini- 
tial value conditions. In particular, the ip-(x,t) com- 
ponent is initially zero, but gradually comes into be- 
ing in the interaction region, as the ip+(x, t) wavepacket 
passes through. The form of the coupling in Eq. (|14[) all 
but assures at least this much. In addition, however — 
and quite unexpectedly — we also find that the ip-(x,t) 
time-evolution undergoes two distinct stages. In the first 
stage, ip-{x,t) stays in place in the interaction region, 
as it grows steadily in magnitude. Once the ip-(x,t) 
integrated probability has grown to roughly the final re- 
flection probability value, the second stage commences, 
in which ip_(x,t) starts dispersing and moving to the 
left, in roughly the same manner as for a free particle 
Gaussian. This two-stage behavior — somewhat reminis- 
cent of fruit ripening on a tree, and then breaking free — 
is clearly evident in Fig. [31 an "animation plot"— of the 
time-dependent p±(x,t) densities. 

Note that the transition from stage 1 to stage 2 occurs 
at a time substantially after the ip+(x, t) peak has passed 
by that of ip_(x,t). The resultant "time delay" is a 
manifestation of the quantum Goos-Hanchen effect i 40 i 41 
The present bipolar approach thus provides a means of 
measuring this effect directly. The Goos-Hanchen time- 
delay effect appears to be closely related to the lack of 
combined flux continuity discussed in Sec. Ill CI Note 
that for this application, the total integrated probabil- 
ity, f_™ (p+ + P— ) dx, is fairly well conserved over time, 
dipping down gradually from the initial value of 1 to a 
minimum value of 0.86, and then increasing to reach the 
final value of 1 again by the end of the propagation, as 
required. 



2. the %p+ "spur, " and its semiclassical interpretation 

One very interesting and unexpected feature is the 
small spur formed on the left side of V>+(x,t) at inter- 
mediate times. This spur remains behind the main ip + 
peak, staying in place above the growing tp_ wavepacket 
during stage 1., and then moving "backwards" with ?/>_ 
during stage 2 (as is more evident in calculations for pa- 
rameters other than those used in Sec. IIII A ip . At a 
certain point in time, either before or after the start of 
stage 2., the spur starts to diminish in size, and eventu- 
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(c) t = 7280 a.u. 
0.3 




(d) t= 11600 a.u. 
0.3 



-10 -5 5 
position, x 



0.2 - 



CZ) 



0.0 



1 I 1 I 1 

t 

II 
II 

1 1 

II 


1 1 


1 1 

1 1 




1 1 

/ V 




I.I.I. 



-10 -5 5 
position, x 



10 



FIG. 2: Wavepacket dynamics for the symmetric Eckart barrier system with m — 2000 a.u. Each plot represents a "snapshot" 
for a specific time, t, as listed (all units are atomic units). Various component wavepacket densities as a function of position 
are indicated as follows: incident /transmitted, p+{x) = \ip+(x)\ 2 , (solid); reflected, p~(x) = \tp-(x)\ 2 , (dotted); total, p(x) = 
|-!/j(a;)| 2 , (dashed). Initial and final densities, e.g. p°(x) and p^(x), are presented in (a) and (d) respectively, in which the 
potential energy is also represented (dot-dashed). In (b), the main peak of ip+(x) has just passed that of i[>-{x), though the 
spur (clearly visible) is forming. The interference is mainly type I, and tp-(x) is in stage 1. In (c), ^-(x) is in stage 2, and the 
interference is type II — caused by the if>+{x) spur, even though it has dissipated almost completely by this point. 



ally dissipates completely. From Fig. [21 the role of the 
spur is clear — i.e., to bring about interference in the left, 
or "reflected," part of ip(x, t), at intermediate times. 

It is well-known that the "reflected" part ofip(x,t) ex- 
hibits nodes and interference to a far greater extent than 
the transmitted part — indeed, one common strategy in 
traditional unipolar QTM calculations is to ignore the re- 
flected part altogether, after a certain point in time>^i To 
some extent, the observed interference in tp{x, t) is due to 



the the incident /transmitted and reflected contributions. 
However, this cannot be the whole of the story, for the 
incident/transmitted contribution [main ip + (x,t) peak] 
often passes by, long before the interference goes away 
completely. Indeed, the so-called "node healing" process 
may continue even after -0_(x,t) has started to move 
away from the interaction region. Eventually though, the 
nodes will be healed completely — leading to a smooth fi- 
nal tjj^(x) reflected branch, and causing the ip+(x, t) spur 
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FIG. 3: Component wavepacket densities, p±(x,t) as a func- 
tion of position, and for a variety of times, t , for the symmetric 
Eckart barrier system with m = 2000 a.u. The upper family 
of curves represent p+(x) at different times, whereas the lower 
family of curves represent p-(x) (magnified by a factor of 4 x ) . 
The motion of the former over time is left-to-right, whereas 
that of the latter is right-to-left. At intermediate t, a station- 
ary p+(x) spur forms and then dissipates in the interaction 
region, corresponding with the simultaneous "birth" and sta- 
tionary growth of p-(x,t). This is "stage 1" of the reflected 
wavepacket dynamics (Sec. IIII A ljl , represented in the figure 
by the right-most four p—(x) curves. A sudden transition to 
stage 2 dynamics is then observed, in which p~(x) moves to 
the left. 



to vanish from existence. 

It would thus appear natural to interpret the above 
as the result of not two, but three separate contribu- 
tions, as considered briefly in Sec. Q] each of which is 
smooth and well-behaved at all times. In particular, 
the ip-(x, t) component is the reflected contribution, the 
main peak of ip + (x,t) is the incident/transmitted con- 
tribution, and the ip+(x,t) spur — non-zero only at inter- 
mediate times — is the third contribution. There are thus 
two distinct types of interference: type I, between inci- 
dent/transmitted and reflected contributions, and type 
II, between spur and reflected contributions. It is quite 
remarkable that the Eq. (114j) decomposition should be 
generally capable of smoothly disentangling such varied 
and complex interference patterns — a very delicate bal- 
ance is evidently required — yet this appears to be the 
case. 

In Fig. [5Jc), for instance, one can discern interference 
in the p(x) plot far to the left of the interaction region, 
where p+(x) itself appears to be vanishingly small. Yet 



this tiny ip + (x) contribution is very significant, for with- 
out it, p(x) would equal p-{x), which by visual inspection 
is clearly false. There is thus a pronounced sensitivity in 
any Eq. ^ bipolar decomposition — owing ultimately to 
the square-root relation between ip and p — which in prac- 
tice, renders it exceedingly difficult to completely disen- 
tangle interference effects. 

Returning to the idea of a tripolar decomposition, we 
comment that further justification for this interpretation 
can be provided using semiclassical arguments. 28 In par- 
ticular, consider an Eckart barrier scattering problem for 
which the initial Gaussian wavepacket is spreading (dis- 
persing) at the initial time to [unlike Eq. This stip- 
ulation is necessary in order that the classical trajectories 
of the ensemble follow different orbits, so that partial 
transmission and reflection can be achieved in a semi- 
classical context. Figure S] represents the time-evolution 
of such a classical trajectory ensemble, or "Lagrangian 
manifold" (LM)^' 24 ' 25 ' 26 i 42 ' 43 ' 44 ' 45 in phase space. 

The semiclassical LM, though initially single- valued, 
develops caustics over time, rendering it multivalued at 
later times. In particular, one caustic (labeled "A" in 
Fig. forms on the left of the barrier peak, essentially 
sweeping through all of the individual trajectory turning 
points, for those classical trajectories with insufficient 
energy to clear the barrier. For a fairly narrow initial 
momentum distribution, this caustic will not move very 
much in x-space, over time. At a given time, those trajec- 
tories that have moved through caustic A constitute the 
reflected branch of the wavepacket, whereas those that 
have not comprise the incident/transmitted wavepacket. 
In any event, caustic A is responsible for type I interfer- 
ence. 

Let us assume that the leading trajectories of the en- 
semble have sufficient energy to clear the barrier — thus 
giving rise to the transmitted branch of the LM at suf- 
ficiently large times. Since the LM must be simply con- 
nected, there must be a second caustic to the left of 
caustic A (labeled "B" ) , which allows the LM to double 
back through the interaction region and connect with the 
transmitted branch, as indicated in the figure. Caustic B 
represents the leading edge of the reflected wavepacket, 
and as such, continually moves to the left. Over time, 
the connecting thread between B and the transmitted 
branch gets pulled apart by the separatrix "like taffy," 
so that the integrated thread probability becomes van- 
ishingly small. The connecting thread itself therefore 
corresponds to the "0+ spur, giving rise to type II in- 
terference, and a tripolar LM representation in the "re- 
flected" part of ip(x,t). The present bipolar approach 
thus achieves classical correspondence in the classical 
limit. In addition — as also observed previously in the 
bipolar treatment of stationary state o 24 i 25 ' 26 i 27 ' 30 — it also 
leads to smooth, classical-like behavior far from the clas- 
sical limit. 
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FIG. 4: Semiclassical Lagrangian manifold (LM) dynamics for the symmetric Eckart barrier system. Each plot represents 
a "snapshot" for a specific time, t. Thin solid curve represents the Eckart potential. One-dimensional array of filled circles 
represents discrete ensemble of classical trajectories, comprising the LM in phase space at each point in time: (a) the initial 
t — to LM is single- valued with respect to x; (b) by t = ti, the LM is still single- valued, but is about to form caustics; (c) at a 
later time t = t2, two caustics have formed, rendering the left side of the LM tripolar; (d) caustic B moves to the left over time. 



3. electron mass 



The general conclusions discussed in Sees. IIII A II 
and IIII A 21 have also been confirmed for a wide variety 
of other parameter value choices for the Eckart barrier 
potential [Eq. ([55)) ] and initial wavepacket [Eq. in- 
cluding deep tunneling applications. In the interest of 
brevity, we forego additional discussion of most of these 
additional calculations. However, in light of the final 
comments in the preceding Sec. IIII A 21 there is one par- 
ticular case that merits further attention — i.e., the quan- 



tum limit in which m — > 0. To this end, we have per- 
formed bipolar wavepacket dynamics calculations for the 
Eckart system using the electron mass m = 1, rather than 
the proton-like m — 2000. Even in this extremely non- 
classical regime, we find that the dynamical characteriza- 
tion of ip±( x i t), as discussed in the preceding subsections 
still holds true. 

In atomic units, the particular parameter values used 
are as follows: Vq — 20; a — 1.0; m = 1; 7 = 1.0; 
x = —7.5; po — ^7.74597. The po value corresponds 
to an incident kinetic energy of 30 a.u., somewhat above 
the barrier peak, but leading to substantial reflection and 
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transmission. Figure [5] shows the resultant wavepacket 
dynamics for il)±(x,t) and tp(x,t) at four representative 
time slices, including t = to = and t = tf = £ max = 
2.5 a.u. In comparison with Sec. IIII A 11 the dynamics 
is of course much faster, and there is substantially less 
interference, as expected. Another difference is that no 
ip+(x, t) spur is evident, so that the interference appears 
to be mainly of the type I variety. In other respects, 
however, the situation is similar to the proton-like mass 
case — in particular, all three conditions of Sec. Ill CI are 
clearly satisfied. This system also exhibits a substantial 
Goos-Hanchen time delay. 



B. Barrier ramp system 

The barrier ramp scattering system has an asymmetric 
potential with a barrier. It serves as a generic reaction 
profile for any direct chemical reaction, and is thus an 
important benchmark system. The potential functional 
form consists of an Eckart barrier added to a hyperbolic 
tangent function, i.e. 

V(x) = V sech(aa;) 2 + f — ~ — j [tanh ((3x) + 1] , 

(36) 

such that the limiting values are V(xl/r) = Vl/r, as ex- 
pected. In atomic units, the particular parameter values 
used for the first barrier ramp calculation are as follows: 
V = .0020; a = p = 2.5; V L = 0; V R = .0008; m = 2000; 
7 = 0.35; x = —7.0; po = 4. These parameters are cho- 
sen to correspond to those in Sec. IIII A 11 except that a 
larger po value is required, in order that Eq. (j!6p still be 
true with an upper limit of p mm rather than (Sec. HID]) . 
The po value chosen corresponds to an incident kinetic 
energy of .004 a.u., which is substantially above the bar- 
rier peak of ^.0024 a.u. 

As per the discussion in Sec. Ill Dl two separate prop- 
agations are performed, using Eq. I]14p for two different 
sets of initial conditions, to a final time i max = 9570 a.u. 
The first propagation is for the left tpL±(x, t) bipolar de- 
composition, corresponding to V = Vl = 0, for which 
= if>°{x) and ip1-( x ) = 0- The second propaga- 
tion is for the right ipR±(x,t) decomposition, for which 
the initial value condition is given by the Eq. ([2"6"[) in- 
tegration, which is computed numerically using standard 
Fourier transform methods. The two sets of solutions are 
then spliced together discontinuously via Eq. (|2"Tj) at the 
dividing point, xjj = 0. 

Figure [S] shows the resultant wavepacket dynamics at 
four representative time slices. The behavior is exactly 
as predicted in Sec. Ill Dt and otherwise comparable to 
that of Sec. IIII A 11 — except that the reflected branch is 
smaller, owing to the larger initial po value. In partic- 
ular, all three conditions of Sec. Ill CI are satisfied, and 
the bipolar components ip±(x,t) are smooth and well- 
behaved — except of course for the discontinuous join at 
x = xdi most evident in the p+(x) plot of Fig.[S]Jb). The 



i/j±(x,t) are also interference-free, though ip(x,t) itself 
exhibits substantial interference at intermediate times, as 
in Sec. IIII A II Though difficult to discern directly from 
the figure, ipL+(x,t) forms a small spur, which serves to 
heal the nodes in the reflected part of ij){x, t), as observed 
for the Eckart barrier, and discussed in Sec. IIII A 21 

As an alternative to the above asymptotic joining pro- 
cedure, Sec. Ill Dl also suggests simply working with the 
ipL±(x, t) solutions throughout all x and t. The resultant 
wavepacket dynamics are still anticipated to be smooth 
and well-behaved everywhere, though Condition 1. will 
no longer be satisfied in the (x > XD,tf) limit. As an 
added benefit, moreover, it should be possible to work 
with initial wavepackets with substantial |t/>°(p)| 2 val- 
ues in the < p < p m ln range, as this contribution 
is problematic only for the ipR±(x,t) components, when 
V R > V L . 

To test this assumption, we have performed ipL±(x,t) 
wavepacket dynamics calculations for a second barrier 
ramp problem, using initial and final wavepacket con- 
ditions identical to Sec. IIII A II — i.e., p — ~ 3.28634 
a.u. and t max = 11600 a.u. All other parameters are 
as described above. The results, presented in Fig. 
are very similar to those of Sec. IIII A 11 except as ex- 
pected. In particular, both components are smooth and 
continuous throughout, but there is a difference between 
p(x > xd) and p+(x > Xd), which is clearly evident at 
large times — even though p-(x > xd) itself is barely vis- 
ible. [Fig. ITJJd)]. This is a further manifestation of the 
square-root sensitivity discussed in Sec. IIII A~2l 

In addition to the above two barrier ramp prob- 
lems, various other parameter choices have been con- 
sidered, including an electron-mass version analogous to 
Sec. IIII A 31 In all cases, the resultant bipolar decompo- 
sitions have been found to be well-behaved at all times. 



C. Two-surface system 

As our final test application, we consider a model mul- 
tisurface system with / = 2 coupled electronic states. 
For this two-surface model, both of the diagonal poten- 
tial energy curves, as well as the off-diagonal coupling 
potentials, are taken to be Eckart barriers, i.e. 

V n (x) = V 22 (x) = Vosech(ax) 2 , (37) 
V 12 (x) = V 2 i(x) = A,sech(ax) 2 . (38) 

In atomic units, the particular parameter values used are 
as follows: V = .0024; D = .00072; a = 2.5; m = 2000; 
7 = 0.35; x = —7.0; po = ^.28634. These parameters 
are chosen to correspond to those used in Sec. IIII A 11 
and also to yield substantial final probability for all four 
components, ip{±(x) and ip2±( x )- 

Note that the above model system conforms to the 
asymptotically symmetric special case (end of Sec. Ill El) . 
Consequently, the choice V° = leads to perfect asymp- 
totic separation for all four components (Condition 3). 
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FIG. 5: Wavepacket dynamics for the symmetric Eckart barrier system with m = 1 a.u. Each plot represents a "snapshot" for 
a specific time, t, as listed (all units are atomic units). Various component wavepacket densities as a function of position are 
indicated; see Fig. [2] caption for explanation. In (b), the ip+(x) peak has just passed that of ip-(x), which is in stage 1. In (c), 
tp-(x) has just entered stage 2, implying a pronounced Goos-Hanchen time delay [also suggested by (d)]. Note the evident lack 
of interference, and of a ij)+(x,t) spur. 



Since the initial wavepacket is incident on surface i = 1, 
the only non-zero initial condition is for ip® , (x) , for which 
Eq. (|34p is used with parameter values given above. The 
four components, ij)i±(x, t) and ip2±(x,t), are then prop- 
agated using Eq. pTj) , to a final time i max = 11600 a.u. 
Only a single propagation is required. 

Figure [5] shows the resultant multisurface wavepacket 
dynamics at four representative time slices. For visual 
clarity, the ip2±(x,t) components are plotted above the 
i/ji±(x, t) — though this is not meant to imply that the V22 
potential is higher in energy. From the figure, it is clear 



that each of the four components is smooth and well- 
behaved at all x and i, and otherwise acts as expected. 
As the incident wave ?/>i+ (x, t) passes through the inter- 
action region, all three scattered components, ipi-(x,t) 
and ip2±(x,t), come into being. All three of these com- 
ponents exhibit the two-stage process of first growing in 
place, and then moving away from the interaction region 
in their respective directions. 

Note that, as in the earlier examples, Tp\+{x,t) devel- 
ops a spur, which provides type II interference (node 
healing) for the "reflected" part of ip±(x,t). The spur 
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FIG. 6: Wavepacket dynamics for the asymmetric barrier ramp system using "spliced" solutions, ij)±(x,t) [Eq. (|2ip ]. Each plot 
represents a "snapshot" for a specific time, t, as listed (all units are atomic units). Various component wavepacket densities as 
a function of position are indicated; see Fig. [2] caption for explanation. In (b) and (c), the vertical bar at x = xd = indicates 
the point where the tj)L±(x,t) and ipn±(x,t) solutions are spliced together; the resultant discontinuity is particular evident in 
the p+(x) plot of (c) (solid curve). 



is particularly prominent for this system [Fig. EJb) and 
(c)]. In contrast, neither of the ip2±{x,t) components 
develops a spur, so that %j)2(x,t) has no type II inter- 
ference. In fact, ip2(x,t) does not appear to exhibit any 
interference at all, even though the ip2±{x, t) components 
do overlap slightly at intermediate times, which could in 
principle lead to type I interference. The reason is that 
the ip2±{x,t) are both growing in place, and therefore 
have the same mean velocity of zero, whereas ipi+(x, t) is 
moving relative to ip\-[x, t). This situation, though eas- 
ily understood within the present bipolar picture, would 



perhaps be difficult to justify in purely semiclassical LM 
terms. 



IV. SUMMARY AND CONCLUSIONS 

This paper represents a turning point in the develop- 
ment of the bipolar QTM methodology. All of the pre- 
vious papers in the series have focused exclusively on 
stationary states — whether bound or scattering states, 
for either continuous or discontinuous potentials. This 
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FIG. 7: Wavepacket dynamics for the asymmetric barrier ramp system using ipL±(x,i) solutions throughout x and t. Each plot 
represents a "snapshot" for a specific time, t, as listed (all units are atomic units). Various component wavepacket densities as 
a function of position are indicated; see Fig. [2] caption for explanation. In comparison with Fig. [6j the component wavepacket 
densities are now continuous everywhere, but do not satisfy perfect asymptotic separation in the (x > xo,tf) limit, as is evident 
in (c) and (d). Also, the reflection probability is greater, due to the smaller value of po- 



was a necessary prerequisite for the present paper, in ad- 
dition to being useful in its own right — i.e., leading to 
the development of extremely efficient and robust algo- 
rithms for ID stationary scattering applications (Ref. [27| 
and Ref. [13). These algorithms have been succesfully 
applied to challenging deep tunneling applications, and 
also to ID reaction path Hamiltonian approximations for 
real chemical reactions such as Cl~ + CH3CI — > CICH3 
+ Cl~. On the other hand, the direct calculation of the 
stationary <\> E states, as functions of position, will never 
be feasible for very large systems, owing to the fact that 



such states are delocalized over a configuration space of 
many dimensions. If such systems are to succumb to 
exact quantum dynamical treatment in a position-space 
representation, it must be via some non-stationary, lo- 
calized wavepacket approach. Moreover, the dynamical 
equations used for the numerical wavepacket propagation 
must be free of any explicit dependence on the 4> E states, 
or even E itself. 

This we have achieved here in Eqs. (fTi]) and (|5Tj) — 
through a sequence of manipulations of the original 
time-evolution equations for stationary scattering states 
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(c) t = 6800 a.u. 
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FIG. 8: Wavepacket dynamics for the symmetric multisurface application of Sec. MI CI Each plot represents a "snapshot" for 
a specific time, t, as listed (all units are atomic units). Various component wavepacket densities as a function of position are 
indicated. Those for surface 1 are grouped together at the bottom of each plot as follows: incident/transmitted, pi+(x) = 
|i/>i + (:r)| 2 , (solid); reflected, p\~{x) = |i/)i_(a;)| 2 , (dotted); total, pi{x) = \ipi(x)\ 2 , (dashed). The corresponding surface 2 
densities are grouped at the top of each plot in similar fashion — e.g., transmitted, p2+(x) — \ip2+(x)\ 2 , via the upper solid 
curve, etc. Surface 2 exhibits little or no interference, in part due to a lack of spurs. In contrast, the %j}\j r {x,i) spur is quite 
pronounced at intermediate times, i.e. in (b) and (c). 



[Eq. ([3])] , and the introduction of reasonable restrictions 
on the allowed wavepacket dynamics (Sec. Ill Cp . The 
intermediate result, Eq. (fTT))) . though not used directly, 
depends on the JWKB quantity (Sec. Ill Bp — thereby af- 
fording a theoretical connection with semiclassical me- 
chanics and the classical limit. The final Eq. (fT4"P is an 
integro-differential equation that depends on the spatial 
integral of t/)±, rather than the integral of IV'iPi or 4'± 
itself. This is quite unusual in the context of quantum 



dynamics — though commonplace, for instance, in soliton 
dynamicsi 46 ' 47 In any case, it is precisely this spatial in- 
tegration that removes all explicit dependence of Eq. (fT5p 
on E and p — the crucial requirement in the derivation of 
Eq. (HI. 

The derivation is more involved than originally antic- 
ipated, owing to the fact that Eq. ([3]) is presented in a 
form that is unsuitable for a wavepacket generalization — 
even though ultimately, Eq. (Q3J) is derived from Eq. , 
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and is thus equivalent to it (at least at asymptotically 
large times). Nevertheless, there are essential differences 
between the two equations, stemming from the fact that 
the coupling contribution in Eq. ([3]) is not the same as 
in Eq. (fT3|) . One important difference is that Eq. (fl4]) 
is independent of the constant potential value, Vq, used 
to define the stationary state trajectories (though the Vq 
value does affect the initial value conditions). Also, the 
TDSE is satisfied at all times (not just asymptotically 
large times). Finally, unlike Eq. ([3]), Eq. (fH)) does not 
conserve total integrated probability [Eq. (fTT)) ]. 

Another important new development is that the bipo- 
lar wavepacket methodology readily lends itself to a tra- 
ditional Bohmian mechanics interpretation — e.g., quan- 
tum trajectories defined via S'± — p± [Eq. 1(15])]. with 
dynamics governed by the quantum potential. The be- 
havior of the bipolar quantum trajectories ensuing from 
Eq. p4[) will be explored in a future publication. In prac- 
tical terms, however, the trajectories are essentially guar- 
anteed to be well-behaved, provided that all three of the 
conditions for well-behaved bipolar wavepacket compo- 
nents, as described in Sec. Ill CI are satisfied for all x 
and t. Condition 3., in particular, is required to circum- 
vent the node problem plaguing the synthetic QTM ap- 
proach, and is therefore essential for application to large 
systems, for which the Eulerian fixed-grid algorithm used 
here (however efficiently implemented) would be unfeasi- 
ble. 

To satisfy all three conditions of Sec. Ill CI for a wide 
variety of applications is decidedly nontrivial, and quite a 
lot to ask of any set of dynamical equations — especially 
with regard to the interference Condition 3., owing to 
the square-root sensitivity discussed in Sec. IIII A 21 In- 
deed, it may be the case that no set of evolution equa- 
tions other than Eq. fTi]) would be capable of achiev- 
ing such a separation for general applications. In fact, 
a great many candidates were explored and rejected, for 
failing to satisfy Condition. 3 (and in many cases, Con- 
ditions 1. and 2. as well)^ The examples of Sec. IIIII 
(and others) nevertheless indicate that the present bipo- 
lar wavepacket approach appears capable of achieving 
this — in both the quantum and classical limits, and even 
when there is a complicated interplay of at least two dif- 
ferent types of interference. These examples also clearly 
demonstrate that virtually all quantum effects that play 
a role in wavepacket dynamics — i.e. dispersion, tunnel- 
ing, and interference — can be easily incorporated. Note 
that for real molecular systems with atomic nuclei heav- 
ier than hydrogen, interference effects can be even more 
pronounced than for any of the examples considered in 
this paper; such cases will be considered in future publi- 
cations. 

The present work thus serves to demonstrate that a 
synthetic bipolar wavepacket QTM approach based on 
Eq. (fT4]) would be widely applicable and numerically 
feasible — using the quantum potential to handle those 
quantum effects that it does best, i.e. dispersion and 
tunneling, and intercomponent coupling to treat inter- 



ference. The numerical algorithm would employ well- 
established unstructured grid techniques, such as local 
least-squares fitting ) 11 ' 16 ! 48 that have already been ap- 
plied successfully to standard unipolar wavepacket QTM 
applications without interference. In the unipolar con- 
text, the chief numerical requirement is the calculation 
of spatial derivatives needed to compute the quantum 
potential or quantum force. The only new numerical re- 
quirement for the bipolar treatment is that of spatial in- 
tegration, to evaluate ^±{x,t) at every time step. This is 
not anticipated to pose severe numerical difficulties, even 
for multidimensional applications, as the generalization 
of Eq. (fT5]) is a line integral (ID), rather than a volume 
integral over all degrees of freedom. 

As a practical matter, it is important that the present 
approach can be generalized for arbitrary asymmetric po- 
tentials, and even multisurface applications, as discussed 
in Secs. lilDl andl llEI respectively. Intriguingly the time- 
evolution equations for the former are unmodified from 
the symmetric case, although the initial and final value 
conditions do depend on the asymptotic potential val- 
ues. From a formal perspective, the fact that multiple 
calculations must be performed in the asymmetric case — 
and then "glued" together discontinuously at the dividing 
point, xd — is perhaps less appealing than the continuous 
V eS approach developed for stationary state applications 
in Ref. [3(1 On the other hand, the discontinuous ap- 
proach is also taken by other standard "dividing surface" 
reactive scattering methods, which essentially posit two 
completely different asymptotic Hamiltonians, one for re- 
actants, and one for products . 34 ' 35 In any event, a gen- 
eralization of the bipolar wavepacket theory using step 
effective potentials V cS (x) = V L + (V R - V L )Q(x - x D ) 
directly (Ref. H^), will be explored in a future paper. 
This will allow direct propagation across the discontinu- 
ity, which in turn, implies that only a single calculation 
need be performed for asymptotically asymmetric sys- 
tems, while still satisfying perfect asymptotic separation. 
Another issue that may be considered in future is the 
generalization for initial wavepackets that do not satisfy 
Eq. (|16p — i.e., that include an initially outgoing contri- 
bution. The nominal difficulty is that such wavepackets 
lead to delocalized \& functions, although this issue may 
be resolved by explicit consideration of the right-incident 
stationary solutions. 

In addition to the various "horizontal" developments, 
to be considered in future publications as described 
above, the next step in the "vertical" or methodologi- 
cal direction — and the subject of the final paper in this 
series — remains the development of a multidimensional 
generalization. In a sense, Eq. (|31[) already provides 
one avenue for multidimensional application — in that the 
multiple surfaces may be regarded as parametrized, dis- 
crete quantum states for all of the "perpendicular" de- 
grees of freedom. This approach could be used, for in- 
stance, to treat rotational degrees of freedom in the con- 
text of a partial wave expansion. In paper VI though, 
we shall address wavepacket dynamics directly on the 
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full-dimensional configuration space. The theoretical 
development required will be seen to be a remarkably 
straightforward generalization of that presented here. In 
particular, only two wavefunction components are still 
required — regardless of system dimensionality — provided 
there is a single reaction path. Moreover, the formal- 
ism can accommodate standard Jacobi-type coordinate 
representations with arbitrary curvilinear reaction paths. 
Paper VI will present results for several multidimensional 
applications including collinear H+H 2 — the first such ex- 
act quantum dynamics calculations ever performed us- 
ing a QTM. Moreover, the resultant multidimensional 
ip± decomposition will be found to satisfy the three all- 
important conditions of Sec. Ill CI 

As a final observation, we note that the bipolar 
wavepacket decomposition as presented in this paper is 
by no means restricted to the Bohmian mechanics con- 
text only. More generally, it ought to be regarded as 
a scattering formalism in its own right, which could in 
principle impact favorably on any of a number of exist- 
ing computational methodologies. Even straightforward 
Eulcrian fixed-grid propagations, for instance, might ben- 
efit from the greatly increased grid spacing and time 
step sizes associated with the component ip±[x, t)'s — 
which are much smoother in general than i(j(x,t) itself, 
particularly in the classical limit. Another idea might 



be to exploit the Gaussian-like properties of the com- 
ponent i/j±(x, i)'s to develop an approximate Gaussian 
evolution scheme, a la Heller Perhaps as few as three 
such "growing Gaussians" would be required — one for the 
incident/transmitted wavepacket, one for the reflected 
wavepacket, and a third for the spur. 
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